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ABSTRACT 

3D modifications to the well-studied 2D flow topology around an embedded planet have the potential to 
resolve long-standing problems in planet formation theory. We present a detailed analysis of the 3D isothermal 
flow field around a 5 Earth-mass planet on a fixed circular orbit, simulated using our multi-GPU hydrodynamics 
code PEnGUIn. We find that, overall, the horseshoe region has a columnar structure extending vertically much 
beyond the Hill sphere of the planet. This columnar structure is only broken for some of the widest horseshoe 
streamlines, along which high altitude fluid descends rapidly into the planet’s Bondi sphere, performs one 
horseshoe turn, and exits the Bondi sphere radially in the midplane. A portion of this flow exits the horseshoe 
region altogether, which we refer to as the “transient” horseshoe flow. The flow continues as it rolls up into a 
pair of up-down symmetric horizontal vortex lines shed into the wake of the planet. This flow, unique to 3D, 
affects both planet accretion and migration. It prevents the planet from sustaining a hydrostatic atmosphere due 
to its intrusion into the Bondi sphere, and leads to a significant corotation torque on the planet, unanticipated 
by 2D analysis. In the reported simulation, starting with a £ ~ r 2/2 radial surface density profile, this torque 
is positive and partially cancels with the negative differential Lindblad torque, resulting in a factor of 3 slower 
planet migration rate. Finally, we report 3D effects can be suppressed by a sufficiently large disk viscosity, 
leading to results similar to 2D. 

Subject headings: hydrodynamics, planet-disk interactions, protoplanetary disks 


1. INTRODUCTION 

The interaction between newly formed planets and the pro¬ 
toplanetary disks, in which they are embedded, is an integral 
part of planet formation theory. In recent years, the field of 
planetary science has made great strides with the discovery 
of thousands of planets and planet candidates by the Kepler 
mission (Borucki et al .J2010 >. The majority of these planets 
have sizes of about 1 to 4 Earth-radii (super-Eart hs), located 
at about 0.1 AU away from their host stars (e.g.|Batalha et al. 
2013| Petigura et al. 20131. This wealth of data has enabled us 


to more thoroughly check the accuracy of our theories against 
observations. 

Despite much effort in the study of planet formation, there 
remain some discrepancies between theory and observation 
that urgently need to be bridged. One example is planet mi¬ 
gration. Through gravitational interaction with the circum- 
stellar disk, a planet can gain or lose angular momentum and 
migrate away or toward its host star. Current planet migration 
theory predicts that Earth-size planets located at about 0.1 AU 
away from their host stars would migrate inward and fall onto 
their host stars within a timescale of only abo ut a few thou¬ 
sand years (Masset & Casol i I20T0j |Paardekooper et al.|2010 
|201 1| >, while the typical lifetime of a protoplanetary disk 
about a few million years. The fast inward migration is also 
known as type I migration. Because we do observe many 
planets at these separations from their host stars, this implies 
our predicted migration rate is orders of magnitude faster than 
can be tolerated by observational data. 

Another example is the accretion of planetary atmospheres. 
The study of how planetary cores of a few Earth-masses (M@) 
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accrete gas from their protoplanetary disks is essential for un¬ 
derstanding how gaseous planets are formed. Current planet 


accretion theory often uses the Bondi radius r B to define the 

extent of an embedded planet’s atmosphere (e.g. 

Pollack et al. 

1996; Ikoma et al. 2000; Rafikov 2006; Lee etal 

tporcr 


GM * 

fa = q—r ~. 
c; 


(1) 


where q is the mass ratio between the planet and the host 
star, G is Newton’s gravitational constant, M* is the host star’s 
mass, and c s is the sound speed of the disk. Treatment of this 
form neglects the effects of the background Keplerian shear 
and the disk’s vertical structure. On the other hand, recent 
results in 3D disk-planet interaction by |Ormel et al.| (j2015bjl 
(hereafter OSK15) found that no gas is’’bound to the planet 
over a timescale exceeding tens of planetary orbital periods. 
If their result applies to Earth-size planetary cores, it might 
prevent their transformation into gaseous giants like Neptune, 
since the core’s envelope may be prevented from cooling, con¬ 
tracting, and accepting more gas from the surrounding flow. 

The shortcomings of current theory may in part be due to 
our lack of understanding about the 3D flow of disk mate¬ 
rial near the gravitational influence of a planet. Both planet 
migration and accretion theories rely on understanding the 
flow topology, yet our current picture of disk-planet interac¬ 
tion is based on a simplified 2D flow, confined to the mid¬ 
plane of the disk, or even a ID, spherically symmetric flow 
in the classical Bondi-Hoyle accretion. Figure [T] shows the 
typical 2D corotating streamlines around a planet located at 
r p = (a, (f> p ). In this picture, there are three separate, well 
defined flow patterns: (i) the disk flow, which can be ap¬ 
proximated as deformed, originally circular orbits around the 
host star, represented by the yellow and green streamlines; (ii) 
corotational or horseshoe flow, which contains the horseshoe- 
and tadpole-shaped orbits of the restricted three-body prob¬ 
lem, represented by the red and blue streamlines; and (iii) 
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streamlines bound to the planet, represented by the magenta 
streamlines. We refer the reader to Ormel[ ( 2013| l for a more 
thorough analysis of the 2D flow topology. The following es¬ 
tablished notions in planet formation theory are tied to each of 
the three flow patterns listed above: (i) density waves are ex¬ 
cited in the disk flow; the angular momentum carried by these 
waves results in the Lindblad torque (also called wave torque) 
acting on the planet; (ii) the corotational flow exchanges an¬ 
gular momentum with the planet as it transfers from outside 
to inside the planet’s orbit at r — a, which results in the coro¬ 
tation torque (also called co-orbital torque or horseshoe drag); 
and (iii) the planet is surrounded by an atmosphere separated 
from the other regions of the disk. 

The 2D picture described above is most applicable to large 
planets, whose Hill radius (or Roche lobe radius): 


At 


-(!)' 


( 2 ) 


where a is the semi-major axis of the planet’s orbit, is larger 
than the local scale height of the disk, ho, in order to jus¬ 
tify the flat disk assumption in the vicinity of the planet. For 
super-Earths, this condition is not satisfied. For example, if 
the planet has 1M®, or q ~ 3 x 1CF 6 , then it has r H ~ 0.01a, 
while Iiq is typically between 0.03 to 0.1a. As a result, one 
cannot assume that the flow topology can be adequately de¬ 
scribed by Figure HI What does the horseshoe flow look like 
above and below the planet’s Hill sphere? How does gas flow 
around the planet and how does it become accreted? These 
are the questions we aim to answer in this paper, specifically 
for super-Earths. 

We use high-resolution global simulations to simultane¬ 
ously resolve the global flow in the horseshoe region and the 
local flow within the planet’s atmosphere. In previous works 
on 3D disk-planet interaction, reduction of a varying d egree of 
migration torques has been found. Linear analysis by (Tanaka 
et al. d2002t) has given a 3D torque weaker by about a factor 
of two than in a 2D disk, in line with the earlier analytical 


estimates of Artymowicz 
lations by D’Angelo et a 


(1993), while fully nonlinear simu- 
P'0031 have found that for planets 


with r H less than 60% of ho, the torque can be up to an order of 


magnitude weaker than 2D predictions; and when D’Angelo 
& Lubow|(j2010| studied the fully nonlinear case for a small 
planet !M ffi ), this time they have found good agreement 


with Tanaka et al. (2002). Overall, it appears that 3D effects 
are more important for the larger than for the smaller plan¬ 
ets, which is counter-intuitive. We will demonstrate in this 
paper that these results can be tied together and explained, by 
studying how the flow topology differs from 2D to 3D. 

In terms of the accretion of a planet’s atmosphere, OSK15’s 
results were obtained from the local 3D simulations of sub- 
Earth-mass planets’ atmospheres, which did not allow them to 
self-consistently connect their atmospheric flow to the global 
disk flow. In this work we will not only establish how a 
planet’s atmosphere fits into the co-orbital flow topology, but 
also probe the interesting regime of planet masses closer to 
the critical core mass limit. 


1.1. Plan of the Paper 

In Section[2]we describe our numerical setup. In Section[3] 
we present the flow topology extracted from our simulation. 
In Section[4]we compute the torque on our planet. In Section 
[5] we demonstrate how disk viscosity can affect the flow pat¬ 
tern and consequently the planetary torque. In Section [6] we 
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Fig. 1.— Streamlines around a planet in 2D, plotted in the corotating frame 
of the planet, which is located at the center of this plot. The background 
Keplerian shear is from bottom to top in the inner disk (r < a), and top to 
bottom in the outer disk (r > a). We call the streamlines approaching the 
planet from the inner disk “inner”, and those approaching from the outer disk 
“outer”. The streamlines are color-coded: yellow and green are the inner 
and outer disk flow; red and blue are the inner and outer horseshoe flow; 
and magenta is the flow that is bound to the planet. The crosses mark the 
“stagnation” points, where the velocity is zero. A third stagnation point exists 
at the location of the planet. This point is irrelevant to our analysis, so we omit 
to label it. The streamlines here are computed from a 2D simulation using 
the same setup and resolution as our 3D one (see Sectionpl, but without the 
vertical dimension, and the planet’s potential is not softened. 


conclude and discuss the implications of our results. 


2. SETUP 

To perform global 3D hydrodynamical simulations of disk- 
planet interaction, we use the Lagrangian, dimensionally- 
split, shock-capturing hydrodynamics code PEnGUIn 
(Piecewise Parabolic Hydro-code Enhanced with Graphics 
Processing Unit Implementation), which has previously been 
used to simulate disk gaps opened by massive planets ( |Fung| 
|et al.|[20T4) >. This code is written in CUDA-C, and runs on 
multiple GPUs (Graphics Processing Units). PEnGUIn uses 
the Lagrangian-re map formulism of the piecewise parabolic 
method (PPM; |Colella & Woodward| 11984} ), which is a 
high-order Godunov scheme, with improvements suggested 
by Blondin & Lufkin| ( |l993| ) included. 

Our simulations are done on a fixed [^] cylindrical grid. It 
spans the full 2 n in the azimuth, 0.7a to 1.3a in the radial 
direction, and 0a to 0.15a in the vertical, where a is the fixed 
planet-star separation. PEnGUIn solves the following Navier- 
Stokes equations in cylindrical coordinate: 


Dr 

d7 


S=- p<v »)• 

--Vp- V0>+ -V-T, 
P P 


(3) 


(4) 


3 While PEnGUIn solves hydrodynamics equations in the Lagrangian 
frame, the solutions are rema pped onto a fixed grid after every timestep. See 


|Colella & Woodward] j 1984} for details. 
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where D/Dr is the Lagrangian derivative, p is the density, 
p the gas pressure, v the velocity, T the Newtonian viscous 
stress tensor which depends on the kinematic viscosity v, and 
<D the gravitational potential of the central star and the planet. 
In the center-of-mass frame. 


0 ) = 


GM* 


Jr 2 + r 2 + 2rr x cos (0 - </> p ) + z 2 
_ qGM, _ 

yjr 2 + r^ - 2rr 2 cos (0 - <p p ) + z 2 + r 2 


(5) 


where r x = qa/( 1 + q) and r 2 = a/( 1 + q) are the star’s and the 
planet’s radial positions, respectively; (f> p - n and <p p are their 
angular positions; and i\ is the smoothing (a.k.a. softening) 
length of the planet’s potential. Recall that M* is the mass of 
the star and q is the mass ratio between the planet and the star. 
We set GM*(1 + q) = 1 and a - 1 so that the planet’s orbital 
frequency Q p = 1 and period P p = 2 tt. For convenience, we 

also denote Vk = y/GM„( 1 + q)/r and Q k = \jGM*( 1 + q)/r 2 
as the Keplerian orbital velocity and frequency. We complete 
our set of equations with an isothermal equation of state: p = 
c\p, where c s is the constant sound speed of the disk. 

Like in Fung et al. (j2014|i, PEnGUIn performs in the planet’s 
co-rotating frame, but uses the conservative form of the angu¬ 
lar momentum equa tion in place of explicitly computing the 
Coriolis force ( |Kley|l998j l. Also, we note that in its most ba¬ 
sic form, Godunov"sdieme suffers inaccuracy in resolving hy¬ 
drostatic equilibrium, because the classical Riemann problem 
is constructed without consideration to the source term, which 
modifies the momentum of the fluid only after the flux of con¬ 
servative quantities has been calculated. This typically results 
in spurious oscillations in the numerical solutions. PPM ac¬ 
counts for this by embedding the source term in the Riemann 


solver, as shown by Equations 2.9 in Colella & Woodward 
(1984). We have tested and verified that PEnGUIn adequately 
maintains a 3D disk in hydrostatic equilibrium in the absence 
of a planet. 


2.1. Planet and Disk Parameters 


which can also be written as rg/liQ. The value q t \, ~ 1 marks 
the division line where a planet becomes sufficiently massive 
to significantly modify the disk structure, which is sometimes 
called the “nonlinear” regime. Because we have q t h = 0.56, 
which is close to unity, we do expect our planet to have some 
weak nonlinear effects. The shortest length scale r s can be 
interpreted as the physical size of the planet. The Earth-Sun 
system, for example, would have r s ~ 0.04rn if the Earth 
were placed at 0.1 AU. Since r s is much smaller than both re 
and r b, we expect it to have little influence on the global flow 
topology, but has more influence on the density profile of the 
planet’s atmosphere. 

The initial density profile of the disk is axisymmetric. It 
follows a power law in the radial direction, and is set to the 
hydrostatic solution of an isothermal gas in the vertical direc¬ 
tion: 

/ r \~ 3 _iL 

p(r,z) = pa I- I e 27,2 , (7) 


where we set po = We impose this initial profile ev¬ 
erywhere in our simulation domain, including the vicinity of 
where the planet will gradually appear. The surface density, 
X, obtained by integrating Equation [7] over z, has the follow¬ 
ing profile: X oc r~ 2 . This surface density profile is inten¬ 
tionally chosen to test a prediction about corotation torque. 
It has been shown in 2D that the corotation torque vanishes 
when the re is zero disk vortensity gradient, or (V Xr)/X — 
constant ( Ward|199T Masset & Ogilvie|2QQ4| Paardekooper| 


et al.||2010p , which, for a Keplerian disk, is precisely when 

X oc 2. Consequently, if we find a significant corotation 
torque in our 3D simulation, it will be a new phenomenon not 
captured by 2D analysis. 

The initial velocity field models a steady disk by setting 
d/dt = 0 in Equation HI and ignoring the potential of the 
planet: v = (v r , rO, 0 ), where 


v r 



dlnp 1\ 

dIn r + 2) ’ 


( 8 ) 


= O k 



/z \ d In p 
r ) d In r ’ 


(9) 


We simulate a planet on a fixed circular orbit embedded 
in a viscous 3D disk. The planet-to-star mass ratio is q = 
1.5 X 1(L 5 , which corresponds to ~ 5M® for a solar-mass star. 
We increase the planet mass gradually at the beginning of our 
simulations, starting from zero to our desired value over 1 P p . 
The relevant length scales in this study are ho, the scale height 
of the disk, r H , the Hill radius within which the gravity of the 
planet dominates, r B , the Bondi radius relevant for accretion, 
and r s , the smoothing length. Since c s is constant, we have 
h = c: s /Q k , which has a radial dependence that goes as h = 
/jo(r/a) 2 , and is equal to ho = 0.03 a at a. Our planet has 
ru = 0.017 a (Equation PV, and rg = qGM t /c 2 = 0.0167a. 
Finally, for r s , unlike in 2D calculations where a non-zero r s 
mainly serves as a way to mimic 3D effects, here we include 
it to avoid singularity. We choose r s = 0.1 r X \, or 0.0017a. Our 
set of parameters therefore gives us the following hierarchy of 
length scales: /; > r n ~ r B > r s . 

Another convenient way to quantify the planet’s mass is the 
dimensionless “thermal mass”: 


The kinematic viscosity of our disk is set to v = 10 fi a 2 O p , 
which corresponds to the Shakura-Sunyaev cr-viscosity coef¬ 
ficient a ~ 0.001 ( |Shakura & Sunyaev||1973| >. This choice 
determines the viscous diffusion timescale across the horse¬ 
shoe region. One can estimate this timescale as t v — w 2 /v, 
where w is the half of the radial width of the horseshoe re¬ 
gion. Viscous diffusion can modify the horseshoe flow if t v 
is shorter than the libration timescale of the widest horse¬ 
shoe orbit, t\ ib = (4a/3w)P p . If we approximate w ~ 2rn, 
then our model gives t v ~ 200 P p which is much longer than 
fiib ~ 40//,. As a result, it is safe to neglect the effects of v is- 
cosity in our analysis of the flow topology. In Section [3~i~| we 
will give an exact measurement of w, and in Section[5]we will 
consider the scenario where t v < f| lh . 

2 . 2 . Boundary Conditions. 

Because we cover the full 2n span, our azimuthal boundary 
condition is periodic. For our radial boundaries, it is impor¬ 
tant to prevent the reflection of density waves launched by 


qth = q 



(6) 


4 Since we do not consider the self-gravity of the disk, this normalization 
has no impact on our results. 
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the planet. To achieve this, we impose wave killing zones in 
r = {0.7a, 0.73a) for the inner boundary and r = {1.27a, 1.3a) 
for the outer. Within these zones, we include an artificial 
damping term: 


ax 

dt 


[X(t = 0) - X(0] 


2c s |r - r km | 
(0.03a) 2 


( 10 ) 


where X stands for all disk variables including p, p, and v. 
/"kin = 0.73a for the inner boundary, and 1.27a for the outer. 
In the vertical direction, as we only simulate the upper half 
the of the disk, so we impose a symmetric condition at the 
disk midplane. At the top we use a reflecting condition to 
ensure that mass does not leak in or out. In practice, the sym¬ 
metric and reflecting conditions are equivalent: we copy disk 
variables across the boundary, and reverse the signs of the z- 
component velocity and force. 


2.3. Resolution 

Since the focus of this work is to study the co-orbital re¬ 
gion, we adopt a nonuniform grid that has enhanced resolution 
close to the planet. Individual grid size is calculated using the 
following formula: 

A.r = A(x - x v ) 2 + Ax mm , (11) 

where x can be any one of r, <p, or z. Ax is the grid size, x p is 
the location of the planet, and Av rnm is the desired grid size 
at the location of the planet. The factor A is found using the 
following relation: 





tan(M4Ax mm ), 


( 12 ) 


where jq, is the distance from the planet’s location to the 
boundary of the grid, which equals to 0.3 in r, n in (f>, or 0.15 
in z ■ N is the total number of grid cells within jq,, which 
equals to 144 in r, 1512 in <f>, or 72 in z. The entire grid 
therefore has 288(r) x 3024(</>) x 72 (z) cells. Ajc m j n defines 
the resolution we desire near the planet, which we choose to 
be re/18 ~ 0.001a. This also corresponds to about 2 cells 
per r s . Figure [2] illustrates how the radial cell size changes 
across the grid. This prescription has the advantage of creat¬ 
ing a high resolution region near the planet that has a roughly 
uniform cell size, without introducing abrupt changes in res¬ 
olution which are prone to generating numerical error. Figure 
[3] demonstrates the level of convergence our resolution has 
achieved with respect to our measurement of the horseshoe 
half-width, w. Reducing our resolution by a factor of 1.2, or 
even 1.44, changes our answer by about 1%. 

At this resolution, we have ~ 63 million cells in total. Using 
3 GTX-Titan graphics cards housed in a single desktop com¬ 
puter, our GPU-accelerated code PEnGUIn runs at a speed of 
0.84 seconds per timestep, or about 140 minutes per planetary 
orbit. 


3. FLOW TOPOLOGY 

We run our simulations for 100/ J p to reach a steady state, 
and then compute the time-averaged density and velo city fi eld 
from 100 to 101 P p to obtain our final results. Section 3.1 will 
describ e the overall size and shape of the horseshoe region; 
Section [372] will focus on the horseshoe turn, where a close 
encounter with the planet occurs, and show how this flow in¬ 


teracts with the planet’s atmosphere; and finally. Section 3.3 
will investigate the aftermath of the close encounter. 


(r-a)/r H 



Fig. 2.— Radial resolution of our grid as described by Lquations | 1 1 | and | 1 2| 
Near the planet’s location, we have ~ 32 cells per ho, or ~ 18 cells per r\\. 


3.1. The Horseshoe Region 


Near the midplane, we expect to find horseshoe orbits, like 
those in Figure [T] If we go to a higher altitude, above the 
planet’s Hill sphere, it becomes unclear what kind of trajec¬ 
tory the gas will take to flow around, or across, the Hill sphere. 
One simple question we can ask is how far up in altitude does 
the horseshoe region extend. 

We use our velocity data to reconstmct the fluid stream¬ 
lines. At a given z, the largest A/j = r-ci that still performs a 
horseshoe turn is defined as w, the horseshoe half-width at that 
z; we will also refer to this streamline as the “widest horseshoe 
orbit”. We measure w at <p = <p p - 1 for the inner flow, and 


<p v + 1 for the outer oncjj We find the two differ by about 
~ 1%, with the inner orbits being the wider one. We consider 
this difference insignificant for the scope of this paper. The w 
we report is the average of the two. 

Figure [3] plots w as a function of z. Remarkably, w remains 
nearly constant in z, even when z reaches 6re or ~ 3/zo- Its 
average value, weighted by disk density, is ~ 1.8r H , and the 
corresponding libration time for the widest horseshoe orbit is 
fiib ~ 43P p . This width is wider than what is expect for plan- 

ets in the linear regime (<7 t h ^ 1). which is w ~ 1.2 a y] aq/ho 
(Masset et al. 2006 1 |Paardekooper & P apaloizou 2009| i. In 
these units, our planet has w — 1.35 a ^Jag/hp, This is con¬ 
sistent with the findings of Masset et al. ( |2006) >, where they 
showed that as q t h increases to unity, thereTs an increase in 
the horseshoe half-width compared to linearly estimated val¬ 
ues. We also perform a series of 2D simulations, with varying 
smoothing lengths, to compare with our result. In Figure [4| 
we see that if one stacks layers of 2D disks, increasing r s to 
mimic the effect of increasing altitude, one will not recover 
the same horseshoe half-width we find in 3D. 

To help visualize these trajectories. Figure [5] plots the 3D 
streamlines of the widest horseshoe orbits. Clearly the horse¬ 
shoe region has a columnar stmcture. It additionally shows 
that halfway through the widest horseshoe turns, as the flow 
crosses the planet’s orbit at r = a, some of it rapidly accel¬ 
erates vertically toward the planet, and completes the turns at 
a significantly lower altitude. In fact, all fluid elements that 
started at z < 2/z 0 are pulled down to z < r H after their horse- 


5 It is also possible to measure w post-horseshoe turn, by tracking the 
streamlines backward in time, bu t we find the flow there too complex for 
a clean measurement. See Section [331 for a discussion on the flow topology 
there. 
























3D Simulations of Low Mass Planets 


5 


z i r H ] 

0 1 2 3 4 5 6 



Fig. 3.— Horseshoe half-width as a function of height above the midplane. 
Z refers to the height of the flow before the turn. The magenta dot-dashed 
curve, blue dot-dash-dashed curve, red dashed curve, and black solid curve 
are results from different simulations, where the resolution is 20% higher 
between each curve. The black solid curve is our choice of resolution. This 
plot shows that our measurement has converged to within 1 %. Also shown for 
comparison are results from 2D simulations with different smoothing lengths, 
at the same resolution as the black solid curve (also see Figure|4j. 
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Fig. 4. — Horseshoe half-width as a function of height above the midplane 
in 3D, and a function of smoothing length in 2D. The black solid curve is the 
same as the black curve in Figure [3] The red dashed curve is from a series 
of 2D simulations. This plot demonstrates that a 3D disk behaves differently 
from a combination of 2D layers. 

shoe turns. 

On one hand, this flow is columnar; there is almost no ver¬ 
tical variation in the planar velocities. This appears akin to 
Taylor-Proudman columns, where dv/dz ~ 0 due to a domi¬ 
nating Coriolis force. On the other hand, the vertical flow di¬ 
rectly above (and below) the planet implies a rapid v z variation 
in z, and a more complex, non-columnar flow structure can be 
seen immediately after the horseshoe turn (Figure |5j; these 
flow patterns do not follow Taylor-Proudman theorem. In the 
following section, we will show analytically that a columnar 
flow structure is expected of a Keplerian disk; at the same 
time, we will account for the non-columnar features we ob¬ 
serve. 


3.1.1. Columnar Flow in a Keplerian Disk 


We rewrite Equation[4]in a rotating frame (such as the rotat¬ 
ing frame of a planet), where u is the velocity in the rotating 
frame and D = (0,0, Q p ) is the frame’s rotation frequency. 


-— \-(u- V)M+2r2xtt + Jlx(J7xr) = -VO—V/J+-V-T, 
dt p p 

(13) 

where the terms on the left-hand side (LHS) are the time- 
dependent term, the advection term, the Coriolis force term, 
and the centrifugal force term; on the right-hand side (RHS) 
are gravity represented by potential O, pressure gradient 
force, and the viscous stress. For a barotropic, steady flow 
with a large Reynolds number (Re ~ \u\Iiq/v » 1), we can 
transform Equation 13 into the vorticity equation by taking 
the curl of both sides, 


(u ■ V)tu = (u> ■ V)u - cu(V ■ u ), (14) 

where uj — V x u + 2 Cl is the total vorticity. The second term 
on the RHS contains V • u, which describes the compress¬ 
ibility of the fluid. We can eliminate this term by combining 
Equation [l4| with Equation[3]to obtain: 

(«■ V)£ = (£- V)u. (15) 


where 


u) Vxu + 2 Cl 

P P 


(16) 


is the vortensity (or p otential vorticity) of the fluid. On the 
LHS of Equation [T5fwe have the advection of £, and on the 
right is the vortextilting term. Consider an incompressible 
flow (p = constant), and the case of a small perturbation in 
vorticity ([2f2| » |V x u\). This gives £ ~ 2 fl/p = con¬ 
stant. Additionally, because fl is non-zero only in the z direc¬ 
tion, the RHS also simplifies, and Equation [15] is reduced to 
du/dz = 0, which is the classical Taylor-Proudman theorem 
stating that there is no vertical variation in the flow. 

For a planet’s horseshoe region, we can apply a similar 
analysis, but with relaxed assumptions. First, instead of the 
incompressibility assumption, we allow the fluid to be com¬ 
pressible, but restrict it to a subsonic flow. Quantitatively this 
means the shortest length scale over which p is allowed to 
vary is h. The local Keplerian shear is supersonic far away 
from the planet, so this assumption also restricts us to a ra¬ 
dial range of r ~ a + h; however, in practice, the Keplerian 
shear does not usually generate shocks in the disk, so as long 
as p is smooth, our analysis can apply to a larger radial range. 
Second, we note that a Keplerian disk does not satisfy the as¬ 
sumption |2fZ| » |V x u\, since |V x u\ — 3f2k/2; instead, 
we assume the vorticity of the disk is mainly in the z direc¬ 
tion, such that oj z » to r , oj t! ,. By our assumptions, the LHS 
of Equation fl5l has a magnitude less than j £|. Therefore in 
component term, Equation[l5]can be estimated as: 
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Since the flow is subsonic, we can further approximate j > 


'ttI , 1 £?|;note that 
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Fig. 5. — Streamlines of the widest horseshoe orbits. The left panel shows the inner flow, and the right shows the outer one. Note that 1) the flow has a columnar 
structure along the horseshoe turn; 2) most streamlines go through a sharp drop in a ltitude half-way through their turns, being drawn vertically to the planet; 3) a 
more complex flow structure is seen near the midplane after the turn; see FigureQojfor a close-up picture of the streamlines there. 


an order of unity to our approximation. Finally, rearranging 
Equation [XT] to [19] our order-of-magnitude analysis gives: 


dit. 


dz 

du A 


£ t 12 
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c s 

< — 
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U) z 
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( 20 ) 

( 21 ) 

( 22 ) 


This demonstrates that in the planet’s co-orbital region, the 
variation of v r and in the z direction is suppressed by a fac¬ 
tor of max on the other hand, the vertical motion 

VKIKI/ 

of the gas is allowed to vary over a length scale as short as h. 

The physical interpretation of this result can be found in 
Equation 15 which states that a vortex line can only be tilted 


as much as advection can carry. If there is little planar vortic- 
ity to begin with, the advection of vorticity will not be strong 
enough to tilt the vortex lines to the in the r and (f> directions, 
and as a result the flow must remain columnar. v z , on the 
other hand, is aligned with the vortex lines, so a vertical ac¬ 
celeration will only lead to the compression or stretching of 
the vortex lines, which is allowed through the compressibility 
of the fluid. 

We now verify whether the assumption u> z » a> r , holds 
in the flow around an embedded planet. In Figure [6] we plot 
the ratio /: 


/ = 


r 

Jo 




Kl 


dj 


fpdz 


(23) 


We find / <s 1 in two regions: one where r > a and <p > <p v , 
corresponding to the starting half of the outer horseshoe or¬ 
bits before crossing the planet’s orbit; and another one where 



Fig. 6. — De nsity -weighted vertical average of the planar-to-z vorticity ratio 
(see Equation |23) , plotted as a function of r and 0. Note that / «: 1 in 
most regions, except for two streams corresponding to the finishing half of 
the widest inner and outer horseshoe turns. 


r < a and (f> < 0 p , which is the starting half of the inner horse¬ 
shoe flow. This is consistent with where columnar structure is 
found (see Figure [5jl. In the regions corresponding to the fin¬ 
ishing half of the widest inner and outer horseshoe turns, we 
find / » 1, indicating the planar vorticity overtakes vortic¬ 
ity in z. This corresponds to the more complex flow structure 
after the horseshoe flow crosses r = a near the midplane (see 
Figure|5]l. In Section|3.3|we will further investigate this aspect 
of the flow. 
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3.2. Flow in the Planet’s Bondi Sphere 


In Section 3.1 we identify a rapid vertical motion in the 
horseshoe flow halfway through the horseshoe turn. This is a 
flow moving toward the planet from directly above (and be¬ 
low) it. Here we investigate how this vertical flow affects the 
planet’s atmosphere. If one assumes the planet’s atmosphere 
has an isothermal hydrostatic structure, it satisfies: 


cb/ 

d" 


d<D 

d" 


(24) 


where // is the enthalpy of the fluid, defined by d /7 = dp/p. 
This gives a vertical density profile at the planet’s location: 


Pstatic(z) = Po exp 


H1 


2 h 2 


sjz 2 + rt) 


(25) 


Figure JTJplots p sta tic together with the density profile we find 
from simulation. We find that near the planet, there is a large 
discrepancy between the hydrostatic solution and our simu¬ 
lation result. Within r# of the planet (the Bondi sphere), the 
density can be an order of magnitude less than p sta tic- This is 
because the gas is not at rest. Figure [ 8 ] shows the streamlines 
in the midplancj^jof the disk, with a color code same as Fig¬ 
ure [I] except for the magenta lines (see caption). Comparing 
to the 2D streamlines in Figure [T| the magenta lines, which 
represents the flow of material from within the Bondi sphere, 
have a qualitatively different behavior. In 2D, the atmosphere 
has closed stream lines bound to the planet; in 3D, there is 
no static atmosphere. Rather, there is a mass inflow near the 
planet’s poles, and a comparable outflow in the equator (Fig¬ 
ure [9}. 

So, similar to the conclusion reached by OSK15 in their 
study for small planets (q± = 0 . 01 ), we find that the flow 
within the planet’s Bondi sphere is not static, but instead cir¬ 
culates with the disk. Moreover, it is clear from Figure[5]that 
the flow in and out of the planet’s Bondi sphere is a part of the 
horseshoe flow, and that the two outflowing streams in Figure 

t are simply the continuation of the inner and outer horse- 
oe turns that has been pulled down from above the mid¬ 
plane. In fact, every magenta line in Figure [ 8 ] can be traced 
back to a horseshoe orbit that originated from an altitude of 
about 0.5 ~ 1 ho. We call this the transient horseshoe flow. 


which we will discuss in depth in Section 3.3 Then, where is 
the planet’s atmosphere? 

There is one region where we do find 3D streamlines that 
do not leave the vicinity of the planet, which is a small sphere 
within 1,5r s , or ~ 0.15 /b, of the planet. Recall that r s is 
equivalent to the planet’s physical size. This means we are 
finding that the planet’s atmosphere is not much larger than 
the pre-defined planet radius. However, it should be noted 
that our simulation grid resolves this region by merely 3 cells, 
so the flow there is not numerically accurate, creating sub¬ 
stantial numerical viscosity that lowers the gas density. The 
explicit viscous force in this region also becomes large as it 
scales as r ~ 2 , adding to the already substantial numerical vis¬ 
cosity. In reality, viscosity on this length scale should be much 
weaker than the disk viscosity, since the typical disk eddy size 
is ~ //o, much larger than r s . With a more accurate and real¬ 
istic treatment, the amount of bound gas may be larger than 

6 The midplane is the bottom boundary of our simulation grid, so these 
streamlines approximate the midplane by following the velocity field of the 
bottommost cells, which are centered on z « 0.0005a, with the z-component 
velocity set to zero. 


what we measure. Other than this region, essentially all of the 
gas within planet’s Bondi sphere are part of a more elaborate 
horseshoe flow. 

Our results share similarities with Tanigawa et ah] ( j2012[) , 
who performed local, isothermal, inviscid simula tions with 
a massive planet that has r H = /z (l (r/ lh = 3), and Ayliffe & 
Bate |(|2012jl, who performed smoothed particle hydrodynam¬ 
ics simulations with 15-33 M ffi planets, taking into account 


self-gravity and r adiation. Figure 5 of |Tanigawa et alT ( | 2012 | l 
and Figure 12 of Ayliffe & Bate| ( |2012| l, each showing the 
mass flux across a sphere of radius 0.3 r H around the pla net, 
can be compared to our Figure |7] [Ayliffe & Bate] (|2012j), in 
their “post-collapse” regime after mass accretion "has slowed 
down, find that influx is concentrated in the vertical direction, 
while outflux is only in the midplane. This is in agreement 
with our findings. For the more massive planet that Tanigawa 
et al. ( 2012 [ > simulated, they find both the influx and outflux of 


mass are more concentrated along the midplane, even though 
their influx is still noticeably offset to a higher latitude. This 
is consistent with our expectation that vertical motion plays 
a lesser role in the horseshoe flow for planets with ru > //q. 
Additionally, Tanigawa et al. ( |2012[ ) reported that, similar to 
OSK15’s results and our simulation, the gas is unbound at 
distances larger than ~ 0.1 r H from the planet. This may be 
indicating that the unbound nature of the gas in the Hill (for 
larger planets with r b > ru) or Bondi (for smaller ones) sphere 
is irrespective of planet mass. 

OSK15 introduced the concept of replenishment timescale, 
Replenish, which measures the total amount of mass within the 
Bondi sphere, Mbs, divided by the influx of mass into it, M ln . 
In our simulation, we have Replenish ~ O p 1 • Since we do not in¬ 
clude any sink cells to treat planet accretion, we expect the net 
mass flux across the Bondi sphere to be zero in steady state. 
We find that M ]n is balanced by M out , the outflux of mass, to 
within 1%, with M ln being the larger one. The Bondi sphere 
is therefore still slowly accumulating mass after 100P p , but 
on a timescale much longer than f re pienish- so it is effectively 
in a steady state. If we lower resolution by factors of 1.2 to 
rH/Av m i n = 15 and 12.5, the flow pattern we see in Figures 
[8] and [9] remains qualitatively the same, but M B s decreases by 
2% and 5%, while M- m increases by 10% and 30%, respec¬ 
tively. This indicates that, not surprisingly, pressure support 
in the Bondi sphere is better established at a higher resolution, 
and our result has converged to within a level of ~ 10%. 

Our measurement of /“replenish is much smaller than the 
Replenish ~ lOOOp 1 reported by OSK15 for their < 7 * = 0.01 
planet. The main difference between our result and theirs, 
is that the vertical motion near the planet is much faster in 
our case, which leads to more kinetic support, and lower den¬ 
sity near the planet comparing to p stat ; c . Overall this gives a 
higher Mj n , smaller Mbs, and shorter f rep ienish- For compar¬ 
ison, OSK15 found a density profile that nearly exactly fol¬ 
lows the hydrostatic solution (see their Figure 4), which in- 
dicates a much s lower motion within the Bondi sphere, while 
D’Angelo et al. ( 2003| ) reported supersonic vertical motion 
above the planet for planet masses > 20M®. This suggests 
to us that the Bondi sphere is increasingly more kinetic sup¬ 
ported as planet mass increases. 

3.3. Transient Horseshoe Flow and Wake Vortices 

By now, it is evident that in 3D, there exists a significant 
asymmetry in the flow pattern around the planet before and 
after the horseshoe turn. There is a new flow, which we call 
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distance from planet [ r g ] 



Fig. 7.— Gas density as a function of distance from the planet. Black solid 
curve plots the vertical density profile; red dash-dot-dotted and orange dashed 
curves are both azimuthal profiles in the midplane, but in the increasing and 
decreasing direction of 0 respectively; similarly, blue dash-dot-dotted and 
green dashed curves plot the radial profiles in the midplane, and are in the 
increasing and dec reasing direction of r. Black dashed curve is calculated 
with Equation 1 25 1 Note the large discrepancy between the black solid and 
black dashed curves. Comparing the black solid curve to the four profiles in 
the midplane, we see that the density structure near the planet is flattened by 
a factor of about 2/3. 



Fig. 8.— Streamlines in the disk midplane. Compare with Figure [11 for 
differences between 2D and 3D flow. Yellow, red, green, and blue streamlines 
are assigned in the same manner as Figure[l] Unlike Figure]!] magenta lines 
are outflows away from the planet, pulled down from initially nigher altitudes. 
They reach as close as 1.5r s from the planet and are unbound. 


the transient horseshoe flow, where fluid at high altitude is 
pulled down toward the planet, enters its Bondi sphere, and 
exits radially in the midplane. These are shown by the ma¬ 
genta lines in Figure [8] for the midplane. 

The word “transient” refers to the fact that this flow only 
performs the horseshoe turn once, even in steady state. Be¬ 
cause of the fall in altitude, gravitational potential energy is 
converted to kinetic, and the flow gains radial speed as it com- 



Fig. 9.— Mass flux across the surface of a sphere centered on the planet. The 
sphere has a radius of 0.5 i'b . Blue and green indicate influx; red and yellow 
are outflux. The speed of the downward flow is about 0.7c s in this plot, 
while the two radial outward flows in the midplane (one not visible from this 
viewing angle) each has a speed of ~ 0.2c s , as is explained in Appendix [A] 
Match this figure with Figure[S]for a more complete view of the flow topology 
near the midplane. 


pletes its horseshoe turn. The transient flow is the portion of 
the horseshoe flow that has gathered enough speed to over¬ 
shoot radially and exit the horseshoe region. We analyze the 
radial outflow in Appendix[A]and show that, due to the conser¬ 
vation of Bernoulli’s constant, the radial outflow at \r-a\ = re 
has a speed of \u T \ ~ 0.6c s , while we measure 0.2c s to 0.4c s in 
our simulation. Appendix [A] also demonstrates that this out¬ 
flow speed is independent or planet mass in the limit r/n, » 1, 
but is slower for lower mass planets. 

Since material in the horseshoe region is being lost to the 
transient flow, it must be replenished. This comes from the 
disk flow lying outside of the horseshoe region; as the tran¬ 
sient flow completes its horseshoe turn, it becomes vertically 
compressed, creating a low pressure region above it, which at¬ 
tracts the high altitude flow right next to the horseshoe region 
to move in (Figure fl0|. This establishes an exchange of ma¬ 
terial between the horseshoe region and the disk , as the red 
and green streamlines in the figure wrap around each other. 

Finally, the vertically compressed transient flow needs to 
decompress as it settles into the disk. This is done through 
a meridional circulation triggered by a “vortex roll-up”, il¬ 
lustrated by Figure TT| The same meridional flow was also 
identified by Morbidelli et al. ( 2014) > in their 3D simulation 
for a Jupiter-mass planet. As the transient flow over-shoots 
the horseshoe region, it gets deflected upward by the mid¬ 
plane disk material, and rolls over to decompress itself, caus¬ 
ing the vortex roll-up. This phenomenon is similar to the be¬ 
havior of the heads of plumes of fluid intruding into station¬ 
ary fluid, which tend to roll at the edges to form mushroom- 
head shapes. However, the closest analogy is the formation of 
wingtip vortices in aerodynamics of finite-span wings, where 
higher pressure air underneath most of the wing length be¬ 
gins to move down (in a downwash), but near the wingtip also 
moves sideways and circulates around the tip to move into 
the low pressure region above the wing. The downwash and 
the associated aerodynamical lift are due to the circulation 
induced around the wing by vortex lines attached to the air¬ 
foil spanwise (Joukovsky theorem states the proportionality 
of lift and circulation). By Kelvin’s circulation theorem, the 
vorticity to, a solenoidal, divergence-free field, cannot simply 
end at the ends of the finite-span airfoil. Instead, the vor¬ 
tex line remains continuous and preserves the circulation; it 
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Fig. 10.— Streamlines at the boundary of the horseshoe region. The red lines are inner horseshoe flow; green are outer disk flow. After a close approach to the 
planet, the red streamlines turn around and descend to the midplane of the disk, sliding underneath the green streamlines. Green lines in higher altitude simply 
enters the horseshoe region, while lower ones are mixed with the red lines. Similarly, but not shown here, this also happens between the inner horseshoe flow and 
the outer disk flow. 



[Cs] 

0.41 


0.31 


0.2 L 



Fig. 11.— Velocity field on a meridional plane at (f> = (f)p-0.5ru/a . The color 
of the arrows indicates the speed. The fastest radial flow speed is ~ 0.4c s in 
this plot. The vortex roll-up occurs between 0.5 ~ lrn away from the planet, 
and about 0.5rn above the midplane. The size of the vortex core is about 
0.1 ~ 0.2rn. 


only changes direction and is shed from the wingtips into the 
wake as two free wingtip vortices separated roughly by one 
wingspan. 

In the protoplanetary disk, we have a close analogy: the gas 
in the transient horseshoe flow after a passage near the planet 
is torqued by its gravity and forms two outflow plumes, one 
directed toward the star and one in the opposite direction. The 
vertical extent of this flow is ~ ra, which plays the role of the 
aerodynamical wingspan of airfoil; two equivalent airfoils of 
this length would be positioned vertically and separated radi¬ 
ally, also by a distance on the order of a few r^. Therefore, 
the planet’s torque on the gas necessarily creates radial plume¬ 
like flows of gas, but also a total of four, alternately turning, 
nearly horizontal vortex lines shed into the disk flow near the 
interface between the horseshoe and disk regions. Figure [T0| 
shows one of these vortices. We call them linear wake vor¬ 


tices, as a reminder that compared with their core diameters, 
they can be very long, as they are carried by the disk flow. 

Furthermore, the equation of vortensity Equation [15] sup¬ 
plemented by viscous dissipation term, which we dropped be¬ 
fore, can be written as 


Dr 


= (£■ V)m + vV 2 £, 


(26) 


where £ is the vortensity given by Equation 16 The evolution 
of vortensity in the core of the vortex can be deduced from 
this equation, after noticing that vortensity and velocity u are 
almost parallel and directed along the vortex line. 

The first term on the RHS of the above equation is thus the 
gradient of longitudinal velocity of material along the vortex 
core multiplied by £ : . As shown in Figure |T0] the roll-up of 
the wake vortices happens near the stagnation region, where 
azimuthal motion with respect to the planet is slow. After the 
vortex core forms, it is carried at an increasing speed into the 
disk flow, therefore the velocity gradient is first strongly pos¬ 
itive, and further from the planet drops to zero, as the flow 
becomes an azimuthally non-accelerated disk flow. While the 
flow accelerates, the full time derivative of vortensity grows, 
strengthening the wake vortex exponentially (in the direction 
of vortex line, spatial derivative of velocity equals D In £/Df). 
Only when the vortex starts freely coasting away from the 
planet, several re behind its roll-up region, the second term 
in the equation representing viscous spreading becomes dom¬ 
inant, diffusing the vortex and weakening the vortensity of its 
core. 

Viscous dissipation of the core has an associated timescale 
equal to the squared diameter dr of the core divided by v, dur¬ 
ing which fluid elements are carried in azimuth by differential 
rotation through a distance A ~ (d 2 /v)Q p re. If d ~ O.lre 
as seen in Figure [IT] ru = 0.017a appropriate for our 5M® 
planet, and v = 10 6 a 2 O p , then the estimated distance is 
A ~ 3ru- This is an estimate of the distance on which viscos¬ 
ity in our simulation will double or e-fold the initial diameter 
by assumption equal to O.lrn. The freely flowing line vortex 
can spread and effectively disappear after a few length scales 
A. This estimate agrees very well with our numerical simula¬ 
tions: vortices are observed to weaken substantially past 10 re 
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behind the planet (Figure [6]). It is intriguing to ask what will 
happen to this vortex if the disk is inviscid, which we will 
investigate in a future paper. 


4. TORQUE ON THE PLANET 


As mentioned in Section [27T] our disk has a flat 2D vorten- 
sity profile, oj/'L = constant, which is predicted by 2D analy¬ 
sis to have a vanishing corotation torque. Since vortensity is a 
conserved quantity along a 2D streamline, a fluid element per¬ 
forming a horseshoe turn will always maintain a density the 
same as its surrounding if the background vortensity profile is 
constant; hence there cannot be a density difference along any 
given r between the inner and outer horseshoe flow, and no 
corotation torque can be generated. As a result, the only re¬ 
maining torque is the differential Lindblad torque, which, by 
the 3D linear calculations by Tanaka et al. ( |2002 1 , is -2.197)], 


where Tq = Soa 4 Q"q ,2 (/io/a) 2 , for an isothermal disk with 
2 oc r 

The above results may not be applicable in 3D for two main 
reasons. First, while the 2D vortensity «/2 can be a constant 
in the disk, the 3D vortensity u)/p depends on z. For instance, 
in our model, oj/p oc r 3/2 in the midplane even though w/2 


is constant. This is because p(z = 0) 


(see Equation 


[7]). Second, and more importantly, vortensity is a conserved 
quantity along a streamline in 2D, but not in 3D, as is evident 
in Equation [15] In fact. Figure [6] already shows that planar 
vorticity is generated after the widest horseshoe turns. This 
section aims to investigate how the corotation torque behaves 
in 3D. 

Before computing the planetary torque in our simulation, 
we first inspect the density structure near the planet. Figure[l2] 
plots the midplane density scaled by the background density 
and with the axisymmetry density about the planet removed: 


A P 




(27) 


where </>' is the azimuth in a polar coordinate centering on the 
planet. The axisymmetry part of the density does not con¬ 
tribute a net torque, so it can be safely removed for clarity. 
One notable problem we can see in Figure [12] is the four¬ 
armed spiral around the planet. It is a numerical artifact due 
to the local Cartesian grid geometry around the planet. This 
problem was also identified by |Ormel et al.| ( |2015a| ). It can 
be reduced by sufficiently resolving r s , which unfortunately 
is not the case for us, since we only resolve r s by about 2 
cells. The four-armed spiral introduces an artificial torque on 
the planet that needs to be removed. We therefore exclude 
the torque contribution from within 0.5/+ of from the planet, 
shown in Figure[l2]as the black circle. 

Now we compute the torque. Figure[l3]plots the torque dis¬ 
tribution dr/dr, which is the amount ortorque on the planet 
by the disk at a given r: 


d T 
~dr 



d~ d (f>. 


(28) 


where T is the net torque on the planet. Figure[l4]plots the net 
torque as a function of time, demonstrating that our measure¬ 
ments have converged with time. Figure |~i5jshows how the net 
torque depends on our choice of the excised region’s radius. 
We find the radius needs to be at least 0.4rB, or about 7 grid 
cells, to fully remove contribution from the non-physical four¬ 
armed spiral, and our choice of 0.5rB safely accomplishes 


that without excessively removing contribution from the disk. 
We further divide the torque into two components: one from 
within the planet’s Bondi sphere (red curve; contribution from 
|r - r p | < /+, represented by the red circle in Figure 121, and 
one from the rest of the disk (blue curve). While the blue 
curve has the characteristic shape of the Lindblad torque dis¬ 
tribution, the red curve is not a well-known feature. If we 
integrate each curve, the red curve gives a torque of +0.507)), 
while the blue one gives -1.27To. The net torque is therefore 
-0.777)). This is significantly weaker than the result from lin¬ 
ear calculation (-2.197))). 

We also perform a 2D simulation to show this result is 
unique in 3D. The 2D setu p is identical to 3D, except we set 
r s = 0.3/io, since in Section |3TT| we find this smoothing length 
produces the best matching horseshoe width. Our 2D torque 
is -2.867)) (Figure [l4]i, compara ble to the 3D value fro m lin¬ 
ear calculation. For comparison, |D’Angelo et al.|(j2003j> also 
found 3D torques are about one order of magnitude weaker 
than 2D when the planet mass is around 10 M®. 

Similar to how the blue curve has two bumps near r - a = 
+/ 20 , and red curve also has two separate bumps in the inner 
and outer disk, but with reversed signs compared to the blue. 
We will refer to this behavior as “torque reversal”. This was 
also seen by |D’Angelo & Lubow| (J2010|> in their Figure 15 
where they simulated planets with masses larger than a few 
M®. This fact, together with Fig ure 1 1 2| provides clues to the 
nature of this torque. In Figure |T2| we marked two stagna¬ 
tion points, where |u| =0. We see that in the outer disk, the 
stagnation point is slightly above </> p , while it is below (b v in 
the inner disk. This is significantly different from the 2D flow 
pattern (see Figure [TJ, where both stagnation points lie much 
closer to <p v . Because of this offset, high density regions are 
created at </> larger (less) than </> p in the outer (inner) disk, thus 
generating a positive (negative) torque near the planet. How¬ 
ever, it remains unclear to us why the net contribution from 
the red curve is positive, i.e., the offset in the outer disk is 
larger than the inner disk, conveniently reducing the net mi¬ 
gration rate. 

We also note that even if we ignore the red curve, the blue 

curve still only contributes a torque of-1.37)), significantly 

weaker than ~ -2.27)) from either linear calculations or 2D 
simulations. This may have to do with the fact that the blue 
curve contains both Lindblad and corotation torque. From 
Figure [13] we can see that a large fraction of the torque dis¬ 
tribution coincides with the horseshoe region. This prevents 
us from distinguishing which part belongs to the corotation 
torque, and which part is the Lindblad torque. However, we 
can measure the corotation torque with a different method, 


separate from Figure 13 
We follow a fluid element’s motion starting at </> = <p v - 1 
for the inner flow, or 0 p + 1 for the outer, until it completes 
its horseshoe turn and returns to its starting azimuth. A/ is 
then the difference in the fluid’s specific angular momentum 
between its start and end points. Combining this with the flow 
rate in the horseshoe region, the corotation torque, 7)- R , is: 


Tcr(z) = 7’cR.i(z) + 7 ’cr, 0 (z) 

= I I p\u$\AI d;' dr 

Ja-W[ *J-z 


J r*a+w o r*z 

I p|i/0|A/ d z 

a J-z 


dr 


<f)=(p p + l 


(29) 

(30) 

(31) 
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where Tqr,{ is the corotation torque due to the inner horse¬ 
shoe flow, Lcr.o is the outer one; w; and w 0 are the horse¬ 
shoe half-widths for the inner and outer flow respectively. 
Furthermore, we can separate the contribution from the tran¬ 
sient horseshoe flow by identifying streamlines that settle out¬ 
side of the horseshoe region. Figure [16] plots the differential 
torque |d7cR/dz| on the left panel and cumulative torque 7 cr 
on the right. Here z refers to starting position of the stream¬ 
line, not where the torque exchange happens. A caveat with 
this method is that our velocity field is time-averaged over just 
IPp, whereas the libration time for these horseshoe orbits are 
much longer. Nonetheless, because Figure p~4] shows the 1 P p 
time-averaged net torque has little fluctuation over a libration 
time, this method should be sufficiently accurate for our pur¬ 
pose. 

We find that the one-sided torques, 7cR.i and 7cr, 0 , each 
has a magnitude of ~ 507o. The transient horseshoe con¬ 
tributes ~ 6% of that torque. 7 CR o is stronger than T C R.i by 
~ 3%, and the net corotation torque is 7cr ~ 1.57)). This 
sufficiently accounts for the difference between our measured 
net torque (-0.777’o) and the expected differential Lindblad 
torque (-2.197’o), which strongly suggests that the corotation 
torque is responsible for both the torque from the red curve 
and the reduction of disk torque in the blue curve. Recalling 
that our disk profile would have zero net corotation torque in 
2D, our result accentuates the difference between 3D and 2D 
torques. The magnitude and sign of the net corotation torque 
is closely related to the asymmetry in the stagnation point off¬ 
sets^ and will be the topic for a future paper. 


D’Angelo & Lubow ( 2010| ) calculated the fully nonlinear 
3D torque on a planet with q ~ 3 x 10 6 and a disk profile 
similar to ours, and found the net torque on the planet to be 
T — -2.297’o, consistent with previous 2D results, implying 
the corotation torque is negligible when 2 oc r - i/2 . Their sim¬ 
ulations differ from ours in two ways. First, our planet is 5 
times more massive, and over 20 times larger in terms of q t h. 
Therefore one can expect our result to deviate somewhat from 
linear calculations. Second, their simulation is in t he re gime 
where t v < f| lh , while ours has t v > t\\\, (see Section [2T| . We 
believe the second point is the main cause for the discrepancy 
between their result and ours. In the next section, we will sim¬ 
ulate a more viscous disk, and check whether we can recover 
their result. 


5. DEPENDENCE ON VISCOSITY 

We increase the viscosity in our model to v = 10 \rQ p 
(a ~ 0.01) and investigate how this will affect the flow pat¬ 
tern and the torque on the planet. We will refer to this case 
as the “viscous” case. Under this setup, the viscous diffusion 
timescale across the horseshoe region t v is ~ 15P p , shorter 
than t \ib ~ 43P p . As a result, most of the gas in the horseshoe 
region can only complete one horseshoe turn before being re¬ 
moved due to the background viscous flow. This prevents the 
flow pattern described in Section [3]from fully setting up. We 
find that in this scenario, the flow has significantly reduced 
vertical variation, and becomes more similar to the 2D case. 

Figure [l7] plots the midplane streamlines similar to Figure 
[T] and [8] but for the viscous case and a 2D simulation with 
r s = 0.3/in (unlike Figure [T] which has r s = 0). In contrast to 
Figure[8] the magenta lines in the viscous case are now inflow¬ 
ing streams entering the Bondi sphere that converges at the 
stagnation point where the planet is located. They correspond 
to the gas attempting to stably orbit the planet, but loses an¬ 
gular momentum (with respect to the planet) too rapidly due 



r-a [ r H ] 


Fig. 12.— The midplane non-axisymmetric density dist ribut ion around 
the planet, scaled by the background density (see Equation Ek The gray 
lines are streamlines in Figure [8| except with the magenta lines omit¬ 
ted. The crosses mark the stagnation points. They are located at {x,y} = 
{-0.47, -0.36} and {0.42,0.42}, in units of rn. This is different from the 2D 
case (Figure [TV where the stagnation points lie close to 0 p . The black circle 
has a radius or 0.5a*b. Because of the non-physical four-armed spiral inside 
the black circle, we exclude this region from our torque calculation. The red 
circ le’s radius is re, corresponding to the sphere where the red curve in Figure 
1 13 | is computed. 


(r-a)/r fi 



Fig. 13.— Torque distribution as a function r. The black solid curve is the 
total torque distribution, and is equal to the sum of the red dashed and blue 
dash-dotted curves. The red curve only in clud es contribution from within a 
sphere of 1 r b around the planet (see Figure |f~2) , while the blue curve includes 
the rest of the disk. The two black dotted lines draw the boundaries of the 
horseshoe region. The two blue bumps at ±ho correspond to the outer and 
inner Lindblad torques; and the two red bumps near the planet are caused by 
the stagnation point offsets seen in Figure [12] 

to both numerical and explicit viscosity. The speed of these 
magenta lines is very slow close to the planet, less than 1% of 
c s , suggesting much of the flow diverges from the midplane 
before ever reaching the planet. As a result, most of the influx 
of mass is still from the vertical direction above the planet. On 
the other hand, the asymmetry about r-a is much reduced 
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Fig. 14.— Net torque on the planet as a function of time. The black curve is 
our 3D torque measurement; red is 2D. This 2D case shares the same setup 
as the 3D one, except for r s = 0.3ho. Both curves are running-time-averages 
over 1 P p . The instantaneous values of the the torques are shown as the gray 
shades around each curve. The vertical dotted lines mark the libration time 
of the horseshoe orbit. The first dotted line is at 1 t \jb = 43 P p , and second 
one is 2 t \jb = 86 P p . 


radius of exclusion [ h Q ] 



Fig. 15.— Net torque on the planet as a function of the radius of the excluded 
sphere centered on t he p lanet, shown as the black solid curve. This plot, 
together with Figure [T2| shows the non-physical four-armed spiral residing 
within ~ 0.4re from the planet contributes a significant amount of torque that 
should be excluded from our calculation. The black dotted line labels the 
rad ius of exclusion we use, 0.5rB, corresponding to the black circle in Figure 

ED 


in the viscous case. As we have shown in Section [33] the 
asymmetry is caused by the vertical flow near the planet, so, 
not surprisingly, we find significantly reduced vertical motion. 
The speed at 0.5 tb above the plan et is 0.1c s , 7 times slower 
than our fiducial case (Section [XT} . 

The more symmetric streamlines also mean the stagnation 
points now lie much closer to 0 p , seen on the left panel of Fig¬ 
ure 17 This reduces the net torque from the horseshoe region. 


Figure 18 plots the torque distribution of the viscous case as 
well as that from Figure [l3]for comparison. It is evident that 
the torque reversal seen before between the blue and red curve 
in Figure [13] has now been erased by the extra viscosity. The 
net torque on this planet is -2.2Tq, while the 2D case on the 
right panel of Figure [IT] has a torque of -2.867o- These are 
in agreement with linear calculations of the differential Lind- 
blad torque in 3D (|Tanaka et al.|2002 1 and 2D (|Paardekooper 


& Papaloizou|2008f, respectively. We therefore conclude that, 
when t v < tut,, the 3D flow field will become similar to 2D, 
and yield a similar torque on the planet as well. 

(D’Angelo & Bodenheimer| ()2013| performed 3D simula¬ 
tions of disk-planet interaction including radiative transfer 
and realistic opacity. They implemented a viscosity of v = 
4 x 10 -6 crQp, and an initial scale heigh j^] of ho ~ 0.06 a at 
the planet’s location. If we assume w ~ 1.2 a \jaq j ho, then 
their model for a 5M® planet is expected to have t v ~ 14/ J p , 
much shorter than fob ~ 70P p . This places their model com¬ 
parable to our viscous case. Comparing the left panel of our 
Figure [17] to their Figure 10, one can see that the flow pat¬ 
terns arelargely similar, and neither of them show significant 
radial outflow at the midplane from within the Bondi sphere. 
This suggests that, in comparison to our fiducial case, the ver¬ 
tical inflow does not penetrate as deep into the Bondi sphere, 
which is consistent with the fact that the inflow speed is also 
much slower. 

Finally, we note that t v > fob is not only a criterion for 
disk viscosity, but also sets a lower limit for the planet mass. 
Namely, it can be rewritten as: 


87 t\ 3 (ho \ I v 


q> 'T) 7 


a 2 Q r 


(32) 


if we approximate w « a Jaq/ho . Together with the fact that 
3D effects are most relevant for embedded planets: ru < ho, 
which can be rewritten as: 


q < 



(33) 


these two criteria bracket the range of planetary mass where 
we expect the 3D flow field to deviate most from 2D. 


6. DISCUSSION AND CONCLUSION 

We present a detailed picture of the 3D flow topology near 
an embedded planet on a fixed circular orbit, extracted from 
3D hydrodynamical simulations of a ~ 5M e planet interact¬ 
ing with a circumstellar disk. Our simulations are run with 
our GPU hydrodynamical code PEnGUIn, on a single desk¬ 
top computer equipped with 3 GTX-Titan graphics cards. We 
found that the 3D modifications to the horseshoe flow have 
a significant influence on both the density structure in the 
planet’s Bondi sphere, and the torque exerted on the planet. 
Below we give a summary of the 3D horseshoe flow: 


(1) At the onset of a horseshoe turn, before the close en¬ 
counter with the planet, the flow is columnar. This re¬ 
sults in a nearly constant horseshoe half-width w in z 
(Figure [3). 

(2) While a fraction of the horseshoe flow continues in 
columnar form after the turn, the widest portion is 
pulled toward the midplane and fall directly on top of 
the planet (Figure [5]> . This flowplummets deep into the 
planet’s Bondi sphere (Figure |9J. 

(3) The release of potential energy from the fall results in 
the flow exiting the Bondi sphere near the midplane at a 

7 There is some ambiguity in the definition of a scale height in their 
work, because of the non-trivial temperature profile in their radiation- 
hydrodynamics treatment. See their Section 4 for details. 
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Fig. 16.— Magnitude of the differential corotation torque on the left panel, and magnitude of the cumulative corotation torque on the right, both as functions 
of height above the midplane. Black sold curve represents contribution from the outer horseshoe flow; red dashed curve corresponds to the inner one. Similarly, 
blue solid and magenta dashed curves are contributions from the outer and inner flows respectively, but are transient flows that exit the horseshoe region after one 
turn. Contributions from inner flows are negative in value. On the left panel, one can see that while the regular horseshoe flow provides the strongest torque near 
the midplane, the transient flow comes from an altitude of ~ ho. On the right panel, it shows that overall the outer horseshoe flow generates a larger torque than 
inner. The sum of all 4 components is 1.57o. 



r-a [ r H ] 



r-a [ r H ] 


F ig. 17.— Midplane streamlines for the 3D viscous case on the left panel, and a 2D case on the right. This 2D case here is the same as the one in Figure 
|14| Comparing to Figure [I] the flow topology in the 3D viscous case is less asymmetric about r = a, and therefore is more similar to the 2D case on the right. 
The stagnation points, labeled as crosses on left, are located at {x,y} = {-0.36,0.08} and {0.36,0.03}, in units of rn- In the 2D case, the large smoothing length 
(r s = 03ho) results in the loss of both stagnation points. Like Figure^ there is a stagnation point at the planet’s location on both the left and right panels, which 
we omit to label. 


speed of order c s (Appendix|A|>. Symmetry of the horse¬ 
shoe streamline about r-a is broken (Figure |8). Con¬ 
sequently the widest horseshoe flow will over-shoot the 
horseshoe region, and exit after just one horseshoe turn. 
We call this the “transient” horseshoe flow. 

(4) As the transient flow pushes into the disk, it is deflected 
by midplane material, resulting in a vortex roll-up (Fig¬ 
ure m At the same time, the loss of material in the 
horseshoe region due to the transient flow is replenished 
by the high altitude flow lying just outside of the region 


Figure [TO] This generates a meridional circulation that 
mixes the flow (Figure [T0|>. 

(5) The meridional circulation (or the ^-direction vortex) is 
eventually killed by disk viscosity, and the flow resets 
to the columnar flow in (1) before the next encounter 
with the planet. 

The flow speed inside the Bondi sphere approaches c s , so 
the gas density there is much less dense than if it has a hydro¬ 
static structure (Figure [7). Nearly all of the gas in the Bondi 
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(r-a)/r B 



Fig. 18.— Torque distribution as a fun ction r. The black solid curve is 
identical to the black curve in Figure [13] The red dashed curve is the torque 
distribution of our viscous case. Note the torque reversal n ear t he planet does 
not exist for the red curve. This is consistent with Figure [17] where we see 
the stagnation points no longer have a large azimuthal offset. 

sphere participates in the horseshoe flow. We found that only 
gas within a distance of ~ r s is bound to the planet. 

We also found that as a part of the asymmetry in the flow 
pattern across r — a, the stagnation points are now offset from 
the azimuth of the planet, (f > v , where the inner point now lies 
below <p p and the outer point above. The flow pattern asymme¬ 
try corresponds to an imbalance in the inner and outer horse¬ 
shoe flow, which generates a net corotation torque of ~ 1.57 n- 
Overall, this results in a net torque of ~ -0.777 q, much re¬ 
duced from -2.197o predicted by linear calculations. 

6 .1. Forming Gaseous Planets 

Following OSK15 to measure the flux of mass in and out of 
the Bondi sphere by f replenish, we found r rep i en ish ~ D p ', which 
is alarmingly short if the Bondi sphere is the planet’s atmo¬ 
sphere. Through streamline analysis, we found that nearly 
all streamlines within the Bondi sphere are not bound to the 
planet. It is then more appropriate to classify the Bondi sphere 
by its flow topology: a part of the transient horseshoe flow. 
However, this leaves us with little atmosphere. If this is uni¬ 
versally true for all planetary cores, then gaseous planets with 
cores of a few M m , cannot be formed. 

There are two major issues with our result in the context of 
gas accretion. First, it should be reminded that we did find 
the gas within 1,5r s of the planet to be bound, but it is only 
resolved by 3 grid cells. Increasing resolution in this region 
will allow us to be more certain about how much gas is truly 
bound to the planet. The resolution required to fully resolve 
the planetary atmosphere is of order the pressure scale height 
on the surface of the planet: /(qGM*) ~ 0.1 r s for our 

setup. This kind of resolution is attainable with local simula¬ 
tions around the planet. 

Second, and more importantly, our simulation is globally 
isothermal, which is unrealistic since it does not take into ac¬ 
count the heating and cooling of the atmosphere. A planet’s 
atmosphere is expected to be heated through the accretion 
of gas and planetesimals, and the timescale for it to cool 
from that heat is much longer than f rep ienish measured from 
our simulation; therefore the planet’s atmosphere should be 
more appropriately described as adiabatic rather than isother¬ 
mal. A more heated and pressurized atmosphere may de¬ 


flect the transient horseshoe flow and resist it from intrud¬ 
ing into the Bondi sphere, allowing more gas to be bound 
to the planet. This may have already been observed in the 
3D radiation-hyd rodynamics simulations by |D’Angelo & B(>] 
denheimer| (|20 1 3), which showed that the planet has a bound 
atmosphere of the size of its Bondi sphere; however, their rel¬ 
atively large disk viscosity (see Section[5]i is expected to slow 
down vertical inflow speed and weaken 3D effects. It remains 
to be seen how large the bound atmosphere is in an inviscid 
flow with realistic radiative transfer. 

6.2. Stopping Type I Migration 

We have shown that our 3D disk produces a net corotation 
torque that is not expected in 2D analysis. As a result, our 
planet, embedded in a 3D disk, migrates 3 times slower than 
if it is driven by the differential Lindblad torque alone. This 
result is subject to a number of uncertainties. 

First, a major limitation to the accuracy of our torque mea¬ 
surement is the poor numerical accuracy near the planet (see 
Figure [T2 ]i. Because of this we have to exclude the r = 0.5r B 
sphere around the planet from torque calculation. Ormel et akj 


showed that numerical convergence can be more eff¬ 
iciently achieved if the grid geometry around the planet is 
polar rather than Cartesian. This is challenging to implement 
in a global simulation, because the grid should, ideally, also 
be polar around the star-planet’s center of mass. Therefore we 
did not attempt it. Nevertheless, we believe our torque mea¬ 
surements do capture the essence of 3D effects, because our 
region of exclusion is sufficiently small that it does not cover 
the stagnation points, and we used an independent method to 
measure the corotation torque which gave a consistent result. 

Second, we have chosen a disk profile that minimizes the 
net corotation torque in order to more easily identify differ¬ 
ences between 2D and 3D. We have seen in Section[4]that the 
one-sided corotation torque has a magnitude that can over¬ 
whelm the Lindblad torque, so an imbalance between the in¬ 
ner and outer horseshoe flow can potentially dominate type I 
migration rate. This can be accomplished by a modification 
as simple as changing the disk density or sound speed profile. 
A future study on how the net 3D corotation torque depends 
on disk parameters will be valuable to understanding planet 
migration. 

Third, we have not considered thermal physics in our 
model, which has been shown to be capable of reversing type 


I migration (e.g. Bitsch & Kley|2011[ 

Bitsch et al. 2014 Lega 

et al.|2014[ Bemtez-Llambay et al. 2015 ). 

Before considering 

ble first step will be 
iiabatic equation of 

the full radiative transfer problem, a possi 
to relax our isothermal condition to an ac 

state. 2D results (Masset & Casoli 2009; Paardekooper et al. 


2010 |> have shown that an additional corotation torque related 


to the enthalpy of the fluid is expected for an adiabatic disk. 
The 3D aspect of this should be studied in more detail. 

Fourth, our planet has a fixed circular orbit, despite the fact 
we are measuring a non-zero net torque on it. From our re¬ 
sults , we can speculate that when the planet migrates inward 
(outward), disk material may be transported from the inner 
(outer) to the outer (inner) disk via the transient horseshoe 
flow, which may further unbalance the inner and outer flow 
and generate a larger net torque on th e planet, leading to type 
III migration ( |Papaloizou et al.|2007 1 . The 3D dynamical in¬ 
teraction between the planet and the disk has the potential 
to modify calculations based on fixed planets by a margin as 
large as the one-sided corotation torque. 

Finally, we note that our planet has a relatively large mass 
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comparing to planets that typically fall in the type I regime, 
which are < 1 M m . This raises the question of how well 
our torque measurement applies to type I migration in gen¬ 
eral. We believe the horseshoe flow asymmetry that causes the 
stagnation point offset is applicable to lower mass planets, be¬ 
cause the flow pattern we presented shares many similarities 
with OSK15’s (compare our Figure[ 8 ]to their Figure 3), who 
simulated the inviscid flow around a planet with q^ = 0 . 01 . 
Therefore, lower mass planets should also experience a re¬ 
duced migration rate like ours. However, it is unclear whether 
the magnitude of the reduction will be similar. This ques¬ 
tion will be answered when global inviscid 3D simulations of 
smaller planets become feasible. 

The ability to simulate smaller planets in 3D is very depen¬ 
dent on computational resources. If we attempt to simulate a 
1 M® planet, then in order to have the same resolution across 
r B as we do in this paper, cell sizes would need to decrease by 
a factor of 5, and the computational time scales as 5 3+1 = 625, 
where the power of 3 comes from the increase in the total 
number of cells, and 1 from the shortened timestep due to 
the Courant limit. This amount of computational resources is 
much beyond what is currently available to us, but in the fu¬ 
ture, advancements in both hardware, such as access to a large 
GPU cluster, and software, such as a more sophisticated grid 
design, may open this pathway for us. 

6.3. Torque and Viscosity 

In Section [5] we suggested to use the condition t v > t \\b to 
identify where 3D effects become important to the flow topol¬ 
ogy and planetary torque. This is also the condition for torque 
saturation in 2D |Ward|[T99T| . Since we do measure a non¬ 
zero net corotation torque, it is important to ask whether this 


torque will saturate. Evidently, Figure [14] does not show any 
sign of torque saturation. 

The process of torque saturation can be described as fol¬ 
lows: because horseshoe orbits are closed in 2D, the entire 
horseshoe region is completely separated from the rest of the 
disk, and therefore has only a finite amount of angular mo¬ 
mentum to exchange with the planet. Consequently, the net 
exchange of angular momentum between the horseshoe re¬ 
gion and the planet in steady state must be zero. The coro¬ 
tation torque can only be unsaturated if fresh supply of an¬ 
gular momentum enters the horseshoe region, which can be 
done through viscous diff usion, hence the torque saturation 
condition. In Section 1331 we showed that there is a constant 
exchange of material between the horseshoe region and the 
rest of the disk due to the existence of the transient horseshoe 
flow. The planet is therefore able to replenish the horseshoe 
flow without help from disk viscosity. In other words, when 
t v > t] it, 3D effects kick in, and the corotation torque is unsat¬ 
urated; when t v < t\ ib, disk viscosity dominates, and it is also 
unsaturated. So corotation torque saturation may be a pure 
2D effect. For a future study, it would be interesting to test a 
range of disk profiles and viscosity in 3D, and find out under 
what circumstance does torque saturation occur. 
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APPENDIX 

RADIAL FLOW SPEED IN THE TRANSIENT HORSESHOE FLOW 


In this appendix we analytically calculate the radial outflow speed at which the transient horseshoe flow exits the horseshoe 
region. Bernoulli’s constant for a given streamline can be written as: 


B = -^r 2 Q 3 + O + \\u\ 2 + q . 


2 p 2 

We divide q into two components: q = ijo + // p , where 770 is the background enthalpy profile that balances the star’s potential: 

d /70 _ d GM, 

dz dj -y/r 2 + z 2 

For x — r - a where x « a, we can rewrite Equation [AT] as: 
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+ ~\ u \ 
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(A3) 


where {x, y, z\ are the local Cartesian coordinates centering on the planet. In Equation A3 770 has been canceled by the z-dependent 
part of the star’s potential. We will also drop rj p because it is expected to have a small contribution to Equation |A3| comparing to 
the planet’s gravitational potential, because as shown in Figure [7] the density, hence pressure, near the planet is much less than 
what is needed to balance the planet’s gravity. 

Consider two points along an inner horseshoe orbit: the first point is just before the fluid is about to perform its horseshoe turn, 
located at pi = {x\,yi,zi}, where xi < 0 , and its velocity can be approximated as u — {u r , 1 , 2 , |a'iQ p , 0 ), where |xiD p is the local 
Keplerian shear; the second point is after the turn, having been pulled down to the midplane, at P2 = U2, >'2, 0 ), where X2 = —x\ 
and yi = >' 2 , and has u = {M r , 2 , |a 2 D p , 0). For convenience we define d 2 - x 2 x +y\- .v 2 +y 2 as the distance from the planet on the 
x — y plane. Equating B at these two points, we have: 


GM, 


„ __— + — 1 — „ 

sjd 2 + z 2 ^ ’ d 

By our results in Section [3TTj it is safe to assume all of the widest horseshoe orbits has the same x\ and y 1 irrespective of their 
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starting z. This allows us to approximate z 2 ~ h^ in Equation |A4| since: 

roo 0 _ _z^_ 

J 0 z 2 p(z = 0)e & dz 


r oo _ 

J 0 p(z = 0)e 2 / 1 “ dz 


= h l 


(A5) 


Finally, combining Equations | A4| and |A5 1 we have: 
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(A6) 


Since Ik^I > |M r ,tl, the horseshoe streamline is asymmetric about x = 0. Equation |A6| should be considered an upper limit for two 
main reasons: first, the streamlines in fact only fall to a fraction of ho instead the midplane; and second, an increase in enthalpy 
as a fluid element moves toward the midplane can compensate for some loss in gravitational poten tial energy. In fact if;; satisfies 
Equation[24] it would leave w r | j = \u T p\, resulting in a symmetric horseshoe flow. Equation |A6| is useful however, for allowing 
us to estimate how this outflow speed scales with planet mass. 

d should scale with the horseshoe half-width w, whi ch can have two possible scalings: w cc a V gja/hp ) for low mass planet s 
( Masset et al.|20 06; Paardekooper & Papaloizou 2009 ), and w cc r H for high mass planets Masset et al.|2006; Peplinski 2008). 


In the low mass limit, if we approximate d 
then we get: 


a Jq(a/ho), and further simplify Equation 


A6 


l«r,2l 


-y 2 V^th' 


cl + Ki| 2 . 


using d <sc h (equivalently, £/, h <sc 1), 


(A7) 


In the high mass limit, we approximate d ~ 2r tt , and have r H » h (equivalently, q± » 1): 


I «r,21 - 



(A8) 


These two limits show that the injection of energy due to the change in z is smaller if planet mass decreases, and asymptotes 
to a constant value for a high mass planet. In our simulation, the radial velocity measured at |x| = ru rang es from 0.2 to 0.4c s , 
depending on the azimuth at which it is measured (see Figures |9] and [TT|>. For comparison. Equation | A8 1 predicts a velocity of 
~ 0.6 c s if | <s c s , which is within an order of unity from our numerical result. 
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